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An imaginary-time nonuniform mesh method for finding eigenvalues and eigenstates of an arbi- 
trary multidimensional Hamiltonian is presented. The main ingredients are (i) a sampling procedure 
optimizing the distribution of grid points and (ii) a diagonalization of a real-valued sparse matrix 
whose eigenvectors are the eigenfunctions evaluated at the selected grid points. The method is 
applied to find the eigenstates of up to five interacting spinless quantum Lennard- Jones particles 
trapped in a ID harmonic potential. Eigenstates of both bosonic and fermionic counterparts are 
obtained, the former exhibiting the phenomenon of fermionization. Finally, the computed excited 
states are used to describe the melting of the Lennard- Jones clusters at finite temperature. 



The multidimensional Schrodinger equation (MDSE) is 
undoubtedly one of the cornerstones of modern physics 
and much attention has been paid to developing effi- 
cient numerical methods for finding its solutions. A 
very rich testing ground for such methods has been pro- 
vided by the development of optical lattices where ul- 
tracold atoms are trapped \V\, and also by the observa- 
tion of new quantum phases at ultracold temperatures 
in finite and homogeneous systems [2l-[7]. Rigorous de- 
scription and explanation of the new physics found in 
these well-controlled experiments require accurate the- 
oretical methods and constitute a formidable challenge 
[S] , the main technical difficulty being the scaling of nu- 
merical algorithms with the number of dimensions D. 
Indeed, standard algorithms for solving differential equa- 
tions, such as the Finite Difference method, scale expo- 
nentially with dimensions [3] , making numerical solutions 
of many-dimensional problems impracticable, if not im- 
possible. Improved methods addressing this difficulty in 
the case of stationary states include the Discrete Vari- 
able Representation (DVR) [HT, collocation method [llj, 
phase-space method based on von Neumann periodic lat- 
tice [T^, variational or diffusion quantum Monte Carlo 
(MC) methods [Hin], Density Functional Theory (DFT) 
[Zl [H] , mean- field or pseudopotential interaction models 
[T51 - frf] . and many others. Some of these methods find 
only the ground state of the time-independent MDSE, 
using different efficient techniques such as the imaginary 
time (IT) propagation [TH| or the Variational Principle 
[19]. Methods for real-time quantum dynamics include 
the Time-Dependent DVR [20,, DFT [21j, mean-field ap- 
proaches [55] , trajectory-based methods such as Bohmian 
dynamics [23 , or time-dependent density matrix renor- 
malization group (t-DMRG) method [24 . The t-DMRG 
has proven to be very efficient for ultracold atoms in op- 
tical traps in one-dimensional geometries, a system ac- 
tively studied due to its interesting dynamical properties 
[551 - E7] . Despite many accomplishments in special cases, 
finding excited states and describing the real-time dy- 
namics governed by a general high-dimensional Hamilto- 



nian remains a difficult computational challenge. 

In this paper we propose a general method, scaling fa- 
vorably with dimensions, which is able to solve the time- 
independent MDSE numerically exactly and simultane- 
ously finds both its ground and excited states. Obviously, 
the proposed IT nonuniform mesh method (ITNUMM) is 
not intended to replace other well established approaches; 
instead we expect it to have a domain of applicability 
where other methods present more technical difficulties, 
such as in finding excited states of many-dimensional sys- 
tems and where efficiency is more important than high ac- 
curacy. To show that ITNUMM achieves these goals, we 
apply it to find the wavefunctions of the first 50 states of 
an ensemble of up to five distinguishable Lennard-Jones 
(LJ) spinless particles trapped in a one-dimensional har- 
monic potential. Once these states are obtained, we 
find, via symmetrization and anti-symmetrization, the 
solutions for the Bose-Einstein and Fermi-Dirac statis- 
tics, respectively, demonstrating explicitly the mecha- 
nism known as fermionization |27| . We also show that 
the computed excited states can be used in a thermal av- 
erage to describe the melting of the LJ clusters at finite 
temperature. 

The derivation of our method starts by rewriting the 
time-dependent MDSE [19, 



ih^^m)) = n\m), 



(1) 



where \ipit)) is the quantum state at time t of the 
D-dimensional system described by Hamiltonian "H, in 
terms of the quantum propagator if(q, q';i — t') :— 
(q|e^*'*^* )'^/''|q') in the position basis |q): 
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Hamiltonian "H := Hq + "Hi is now split into two com- 
ponents: "Ho is any Hamiltonian that includes the ki- 
netic energy operator T and whose matrix elements in 
the q-representation are known, while "Hi = Hi(q) is 
any many-body potential depending only on q. For 



very short time intervals t — t' ~ At, the time evolu- 
tion operator can be split to first order as e"*^*^/'' = 
2Ami/h _,_ (9(Ai2) and one can write 



g-iAt«o/Sg 



7/;(q, t' + Ai) = / dq'Xo(q, q'; Ai)e'''^*^i(i')/'"V(q', t') 

+ 0{At'), (3) 

where KQ{q^, q' ; At) :— (q|e^*'^*^°/''|q') is the propaga- 
tor of Ho, which is assumed to be known explicitly. 
The |q) basis is discretized as 
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where w(qj) is a weight function depending on a partic- 
ular realization of the N states jq^). Indeed, w is defined 
as u;(q) := [Np{q)]^^ , where p(q) is the density distribu- 
tion of the qj. With this discretization, Eq. (pi) becomes 
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Since our main interest is finding the stationary states 
of H, in the following we will assume that (i) '0(q, i) = 
g-ztEn/Sj^^jlq^ where (/9„(q) and _E„ are the nth eigen- 
state and eigenenergy of the Hamiltonian H, and that (ii) 
the evolution is performed in IT (t — >■ —ir). Although 
the density p(q) is arbitrary, below we show that Eq. ((sl) 
simplifies in the IT scheme if this density corresponds to 
the classical Boltzmann distribution of "Hi, namely if 



p(q) = z^] 



,-AT«i(q)/R 
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where Z-^^ = Tre" ^^i''* is a normalization constant 
(called configuration integral) and At jh plays the role of 
the inverse temperature /3. Under these conditions, Eq. 
^h reads 
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By defining vector $„ :— {ifn{(lj)}f=i, whose jth com- 
ponent is the wavefunction evaluated at position q^ , and 
matrix Kjk '■— Ko{qj,qk',—iAT)Z-Ui/N whose elements 
are proportional to the propagator Kq from q^ to q/c , one 
can rewrite Eq. (l7| as a matrix eigenvalue equation 



-ArE^/h^ ^K-^r 



(8) 



In this way, the problem of finding the spectrum and 
eigenfunctions of the original Hamiltonian "H is reduced 
to sampling the classical Boltzmann distribution and di- 
agonalizing the real-valued matrix K evaluated at those 



points. In the special case of T-Lq = T, 'Hi(q) equals the 
classical potential energy, Kq is a free-particle propagator 
in D dimensions [19], and matrix elements Kjk assume 
the Gaussian form 



K,k^ 
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where m is the mass, for simplicity assumed to be the 
same for all degrees of freedom. In correlated sys- 
tems, where sampling the Boltzmann distribution is dif- 
ficult or unfeasible — as in the case of Coulomb interac- 
tion, we propose the splitting 'Hi(q) = Vi(q) + V2(q), 
where Vi(q) is a sum of well-behaved one-body poten- 
tials and V2(q) is the remainder including all correlations. 
Here the sampling is performed with the weight p(q) = 
Z^^^e-^^^i(i)/'' and normalization Zy, = Tre^'^^^i/'"'; 
the matrix to be diagonalized becomes 
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We have found this method to be very efficient in one- 
dimensional problems with several very different poten- 
tials. Although an arbitrary sampling procedure can be 
used, we have employed a quadrature scheme: instead 
of random sampling of p(q) by a MC procedure, the qj 
points are chosen with a deterministic algorithm. The 
motivation for this approach is reducing to a minimum 
the number of vector-elements needed for a given accu- 
racy, and thus reducing the computational cost of the di- 
agonalization of K. Specifically, we first consider a new 
variable u, uniformly distributed in the interval [0, 1], and 
define an equidistant grid Uj = {j — 1/2) /N. The Jaco- 
bian of the transformation from g to m is given by p{q) 
since p{q)dq = du, hence 



u{q) 



dq'p{q') = P{q), 



(11) 



where P{q) is the cumulative distribution function. Next, 
the g-grid is obtained by inverting this equation for all 
values of Uj , and once the g-grid is ready, the evaluation 
and diagonalization of the matrix K is performed with 
standard numerical methods. 

As the first application of ITNUMM, we solved (i) 
the ID harmonic oscillator [19] T-Li{q) = muPq^/2, us- 
ing natural units for energy and position (defined by 
hijj and ^Jh/rnuj, respectively), and (ii) two particles of 
equal mass m interacting via a LJ potential 'Hi{q) = 
Flj(<j) ee e[(re/g)i2 - 2{rjq)% For the latter, we 
used a de Boer quantum delocalization length [28] of 
A — 2^/^h/{rf^^Jrae) — 0.16, corresponding to hypo- 
thetical particles with properties between para-hydrogen 
— where quantum effects dominate — and neon — where 



quantum effects are present but classical behavior dom- 
inates. In the Supplementary Material (SM), we show 
the grid points, eigenvalues, and several eigenstates ob- 
tained with ITNUMM in both cases — we also include a 
notebook executable in the Wolfram Research's Mathe- 
matica software, where the interested reader can explore 
the technical details of the method. 

As expected, we have observed that the imaginary time 
Ar must be small enough to reduce the relative error a 
introduced by the splitting of the time operator — which 
is cr ~ 0(At'^) since the second order term vanishes for 
stationary states — but large enough to avoid reducing 
the Gaussian elements of the K matrix to delta func- 
tions and eventually obtaining a diagonal matrix. The 
latter condition is ensured by requiring 1 3> Kjj = 
Zhi {m/2T:hAT) ' /N, which imposes a lower bound on 
Ar for a given N . In the SM, we explore the dependence 
of the relative error a on Ar for a given number N of 
grid points, and also the dependence of a on N in the 
harmonic oscillator. Remarkably, the relative error can 
be fitted to a{N) ~ 0.18A^~^'^, indicating a significantly 
faster convergence rate than the rate expected for a MC 
scheme [a{N) ^ N''^^^] [3 US]. Regarding the excited 
states, we have found that the error becomes large for 
the states with the highest eigenenergies. Indeed, the 
number of grid points N becomes insufficient to repro- 
duce the characteristic high frequency oscillations of the 
wavefunctions describing highly excited states. Yet, the 
agreement with exact results is very good for the first 150 
states using N = 500 grid points, as shown in Fig. [TJfor 
the first 50 states (the whole spectrum can be checked in 
the SM). 

As a more stringent test, we now apply the method 
to D LJ particles in a one-dimensional harmonic trap. 
Potentials Vi and V2 are defined by 
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^2(q) = E^Lj(kA-g;.|), 
X<fi 



(12) 
(13) 



the de Boer length has the same value as in the example 
above, and uirey/mje = 1/2. The problem is separable 
only for D = 1 ot 2, and so a multidimensional numeri- 
cal method is mandatory for D > 3. In order to reduce 
the number of grid points in the numerical calculation, 
we first solve the problem for distinguishable particles 
and construct a posteriori the eigenstates of indistin- 
guishable particles by symmetrizing or anti-symmetrizing 
the wavefunction for spinless bosons or fermions, respec- 
tively. Thanks to the repulsive nature of the LJ potential 
at short distances we only need to evaluate K in the sub- 
space defined by gi > q2 + a, . . . , qo-i > to + a, where 
a is the core radius of the LJ potential, within which 
the wavefunction is expected to be zero within numeri- 
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Figure 1. (Color online) Left: Energy spectrum for D LJ 
particles in a ID harmonic trap obtained with our method 
(circles). The exact results for D — 1 and 2 are shown as 
red solid lines. Energies are shifted to the minimum of the 
potential min('Hi) — —e{D — l)D/2. Right: One-body den- 
sities (normalized to the number of particles) of the ground 
state for distinguishable (colored lines) and indistinguishable 
particles (black lines). 



cal accuracy (a = 0.63re in our calculations). The grid 
points are sampled from the classical Boltzmann distri- 
bution of the harmonic trap in this subspace, p(q) — 

2;-le-Arrnu^\ci\^/2h ^-^.j^ ^y, = {2tT / ATmUJ^)^/^ Coia) , 

where the normalization constant obeys Cd{Q) = Dl. All 
the two-body interactions, contained in V2(q), are evalu- 
ated in the matrix elements of K. As mentioned above, 
only the low-lying eigenstates are accurate, so we have 
used the Arnoldi algorithm j^S] to obtain the first 50 
eigenstates. We have taken into account that many of 
the matrix elements are close to zero by using standard 
computational techniques for sparse matrices: instead of 
storing the N x N values of the matrix, only elements 
larger than a certain threshold were stored. Parameters 
used in calculations with varying D were 
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500 6709 14 394 36 517 84 690 
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Note the relatively low total number of grid points needed 
to obtain results with reasonable accuracy (a relative er- 
ror of 0.002 for the D = 2 case). Figure [2] shows the 
ground and 19th states for D = 2 and for the three statis- 
tics: distinguishable particles (in the above mentioned 
subspace), bosons, and fermions (in the full space). The 
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Figure 2. (Color online) Wavefunctions </p„(gi,q2) of the 
ground and 19th states for D = 2 LJ particles in a ID har- 
monic trap (see text for details). Lighter (darker) color in- 
dicates positive (negative) values of the wavefunction. Left: 
distinguishable particles in the subspace gi > 92; center: in- 
distinguishable bosons; right: indistinguishable fermions. 




Figure 3. (Color online) One-body densities (normalized to 
the number of particles) for D = 4 LJ particles in a ID har- 
monic trap at three different temperatures: fta;/3 = 2.8 (dotted 
line), 0.7 (dashed line), and 0.4 (solid line). The crystal struc- 
ture disappears with increasing temperature, resulting in an 
unstructured total density as in a fluid. 



spectrum oil-L as a function of D is shown in Fig. [T| (left 
panel). We find the same spectrum for the three cases, 
which is a consequence of the fermionization |27| mecha- 
nism due to the repulsive behavior of the LJ potential at 
short distances. Indeed, the bosonic and fermionic sys- 
tems show the same one-body densities in position space, 
as shown in the right panel of Fig. [T] In all three cases the 
densities show a well-defined structure, forming a quan- 
tum crystal. The displayed one-body densities, defined 
as EHl 



Pn{q\) 
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\<Pn{(i)?X{dq^, 
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were obtained from the nonuniform mesh as follows: 
first, we computed its Fourier transform in a regular 
equidistant grid in momentum (fc) space as 
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p„(fc)=y |^„(q)pe-''=«^I]'^'?' 
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(16) 
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and then Fourier-transformed /0„(fc) back to g^-space us- 
ing standard numerical methods. 

The 50 states obtained in the course of the diagonal- 
ization are sufficient to study the behavior of the system 
at finite temperatures. The (unnormalized) probability 
distribution of the system pp (q) at finite inverse temper- 
ature /3 is defined as the thermal average 



pMq) = ^e-^^"|^„(q) 



(17) 



and the corresponding one-body density Pj3{q) is obtained 
in an analogous fashion as for pure states. Figure |3]shows 
the one-body density for Z? = 4 at three different temper- 
atures, where the lack of structure at the highest temper- 
ature can be understood as the melting of the quantum 
crystal. 

To summarize, we have presented compelling evidence 
that the proposed method achieves the original goals. 
Indeed, (i) the only approximation used is the numeri- 
cal discretization of space and time; (ii) the ITNUMM 
only requires standard methods for sampling from an 
arbitrary probability distribution and for diagonalizing 
real-valued sparse matrices; (iii) both ground and ex- 
cited states are obtained in the course of the diagonal- 
ization; and (iv) due to the nonuniform nature of the 
grid that uses the potential to guide the sampling, the 
complexity of the algorithm is significantly reduced in 
high-dimensional systems. It is worth mentioning that 
the accuracy of ITNUMM can be improved by using a 
splitting of a higher order than in Eq. (Isl . In addition to 
computing thermal averages — as shown here — the large 
set of excited states can be also used for solving real-time 
quantum dynamics in a straightforward fashion. As we 
have not found any a priori limitation to the applicabil- 
ity of the method, other systems described by the MDSE 
will be studied in the future. 

Acknowledgments. The authors thank E. Zambrano, 
M. Wehrle, and M. Sulc for useful discussions. This re- 
search was supported by the Swiss NSF NCCR MUST 
(Molecular Ultrafast Science & Technology) and by the 
EPFL. 



n=l 



* |alberto.hernandodecastro@epfl.ch| 



[i; 
[2: 

[3] 
[4] 

[5] 
[61 
[7] 

J8j 



jiri.vanicek@epfl.ch 

I. Bloch, Nature Physics 1, 23 (2005) | 

P. Kapitza, Nature 141, 74 (1938)] 



[10] 



[18] 



A. J. Leggett, |Phys. Rev. Lett. 29 , 1227 (1972)' 

S. Grebenev, B. Sartakov, J. P. Toennies, and A. F. 

Vilesov, Science 289, 1532 (2000), 

J. R. Anglin an d W. Ke tterle, Nat ure 416, 211 (2002)"] 

S. Balibar, Nat ure 464, 176 (2010)] 

J. M. Mc Mahon, M. A. Morales, C. Pierleon i, and D. M. 

Ceperley, |Rev. Mod. Phys. 84, 1607 (2012)[ 

I. Bloch, J. Dahbard, and W. Zwerger, IRev. Mod. Phys. 



80,885(2008), 
K. Morton and D. Mayers, Numerical Solution of Par- 
tial Differential Equations, An Introduction, 2 ed. (Cam- 
bridge University Press, 2005). 
D. O. Harris, G. G. Engerholm, and W. D. Gwinn, J. 



Chem. Phys. 43, 1515 (1965) 



[11[ R. Kosloff, J. Phys. Chem. 92, 2087 (1988) 

[12[ A. Shimshovitz and D. J. Tannor, Phys. Rev. Lett. 7, 

r~070402 (2012) 

[T3]T. Pollet, |Rep. Prog. Phys. 75 , 094501 (2012)| 

[14| R. G. Parr and W. Yang, Density- Functional Theory of 

Atoms and Molecules, 2 ed. (Oxford University Press, 

1994). 
[15[ F. Dalfovo, A. Lastri, L. Pricaupenko, S. Stringari, and 

J. Treiner, Phys. Rev. B 52, 1193 (199 5) _ 

[16[ L G. Saveiiko, T. C. H. Liew, and L A. Shelykh, |Ph^ 



Rev. Lett. 110, 127402 (2012)| 

J. Rogel-Salaz ar, |Eur. J. Phys' 34, 247 (2013)| 



[19] J. J. Sakurai, Modern Quantum Mechanics , 1 
(Addison- Wesley, 1993). 



G. Billing and S. Adhikari, |Chem. Phys. Lett. 321, 197 

(2000) 



V. S. Popov, |Phys. Atom. Nuclei 68, 686 (2005) 



^0] 

[21] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997| 
L_ (1984) 
[22] M. Abid, C. Huepe, S. Metens, C. No re, C. T. Pham, 

L. S. Tuckerman, and M. E. Brachet, [Fluid Dyn. Res.j 
I 33, 509 (2003) 
[23] R. E. Wyatt and C. J. Trahan, Quantum Dynamics with 

Trajectories: Introduction to Quantum Hydrodynamics, 

1st ed. (Addison-Wesley, 2005). 
[24] G. Alvarez, L. G. G. V. D. da Silva, E. Ponce, and 

E. Dagotto, Phys. Rev. E 84, 056706 (2011) 
[25| J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodg- 

man, S. Langer, L P. McCuUoch, F. Heidrich-Meisner, 

L Bloch, and U. Schneider, [arXiv:1301 .5329 (2013) 
[26[ M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, 

B. Rauer, M. SchreitI, LMazets, D. A. Smith, E. Demler, 

and J. Schmiedmayer, [Science 337, 1318 (2012) 
[27] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Foiling, 

L Cirac, G. V. Shlyapnikov, T. W. Hansch, and L Bloch, 

'Nature 429, 277 (2004) , 
[28[ J. Deckman and V. A. Mandelshtam, |J. Phys. Chem. A 

26, 7394 (2009) 
[29[ W. E. Arnoldi, Q. Appl. Math 9, 17 (1951). 
[30[ E. Lipparini, Modern Many-Particle Physics: Atomic 

Gases, Nano structures and Quantum Liquids (World Sci- 
entific, 2008). 



